Modelo SUR ou MQO? Quais as vantagens e desvantagens?
Modelo SUR
Os modelos de regressões lineares aparentemente não correlacionados ou modelo SUR (Seemingly Unrelated Regression) é um modelo multivariado, isto é, busca explicar múltiplas variáveis dependentes. Além disso, pode ser definido através de \(g\) modelos lineares indexados por \(i\), dados por:
\(\bs y_\bu\) é um vetor de tamanho \(ng\) com os vetores \(\bs y_1, \dots, \bs y_g\) empilhados;
\(\bs X_\bu\) é uma matriz diagonal com blocos \(\bs X_1, \dots, \bs X_g\);
\(\bs\beta_\bu\) contém os vetores \(\bs\beta_1, \dots, \bs\beta_g\) empilhados;
\(\bs u_\bu\) é um vetor de tamanho \(ng\) com os vetores \(\bs u_1, \dots, \bs u_g\) empilhados.
Suposições
\(\E(u_i \mid \bs X_i) = 0\), \(\forall i\), isto é, fracamente exógenos (então, por que não utilizar MQO? 🤔);
\(\E(\bs u_i \bs u_i') = \sigma_{ii} I_n\) e \(\E(u_{ti}u_{tj}) = \sigma_{ij}\), ou seja, a variância do erro é constante em cada equação, mas pode variar entre equações;
\(\E(u_{ti}u_{sj}) = 0\), \(t \neq s\) e \(i \neq j\).
onde \(\sigma_{ij}\) é o elemento \((ij)\) da matriz \(\Sigma_{g \times g} > 0\), denominada de matriz de covariância contemporânea.
Portanto, fica definido o modelo SUR por meio das \(g\) equações em \(\eqref{eq: modelo sur}\) e das suposições apresentadas.
Suposições
Correlação contemporânea?
A correlação contemporânea é quando existe correlação entre os erros dos diferentes modelos para um mesmo indivíduo. Seu surgimento é devido à omissão de variáveis, já que as informações presentes nas preditoras que explicam a variável resposta (mas não entram no modelo) passam a fazer parte do erro.
Vantagens e Desvantagens
Vantagem: Com a adaptação feita, o modelo passa a considerar a correlação presente no erro, aumentando sua precisão.
Desvantagem: Quando o número de indivíduos é muito grande, o número de parâmetros a ser estimado também será, podendo levar à estimativas pouco fidedignas do modelo.
Estimação
Produto de Kronecker
Seja \(A\) uma matriz \(p \times q\) e \(B\) uma matriz \(r \times s\). Então, o produto de Kronecker, denotado por \(\otimes\), entre \(A\) e \(B\) é dado por: \[
A \otimes B = \begin{bmatrix}
a_{11} B & \dots & a_{1q}B\\
\vdots & \ddots & \vdots\\
a_{p1} B & \dots & a_{pq}B\\
\end{bmatrix}
\]
Se \(\sigma_{ij} = 0\), \(\forall i \neq j\), então não há ganhos em aplicar outro método que não seja MQO;
Se \(\bs X_i = \bs X_j\), \(\forall i,j\), então MQO e MQG são iguais;
Quanto maior a correlação entre \(u_{ti}\) e \(u_{tj}\), então maior a eficiência de MQG sobre MQO;
Quanto menor for a correlação entre \(\bs X_i\) e \(\bs X_j\), \(i \neq j\) maior serão os ganhos em utilizar MQG do que MQO.
MQGF
Contudo, no geral, \(\Sigma\) não é conhecida. Então, por MQGF:
\[\hat\Sigma = \frac{1}{n} \bm{\hat U}'\bm{\hat U}\] onde \(\bm{\hat U} = \begin{bmatrix} \bs u_1 & \cdots & \bs u_g \end{bmatrix}\), isto é, uma matriz \(n \times g\) cujas colunas são os vetores de resíduo \(\bs{\hat u}_i\) da \(i\)-ésima equação, obtidos por MQO. Além disso, utilizar \(n\) produz um estimador viesado, então, se \(k_i = k_j = k\), \(\forall i\neq j\), pode ser utilizado \(n-k\) para obter um estimador não viesado para os elementos da diagonal de \(\hat\Sigma\).
Mínguez, R., López, F.A., & Mur, J. (2022). spsur: An R Package for Dealing with Spatial Seemingly Unrelated Regression Models. Journal of Statistical Software, 104(11), 1–43. doi: 10.18637/jss.v104.i11.
Contexto
O SUR espacial é uma técninca que incorpora simultaneamente efeitos espaciais e erros correlacionados entre as equações. Sua primeira menção ocorre em Anselin (\(1998\)): “Uma equação para cada período de tempo, a qual é estimada para uma seção transversão de unidades espaciais …”.
O artigo utilizado para a análise disponibiliza o pacote spsur no R, o qual foi desenvolvido para considerar modelos de multiequações com efeitos espaciais e correlação residual.
Modelo
O modelo geral SUR com dependência temporal SUR-GNM (general nesting model) é dado por:
com \(\bs u_t = \lambda_t \bs W_t^u \bs u_t + \bs\varepsilon_t\), \(\E(\bs\varepsilon_t) = 0\) e \(\E(\bs\varepsilon_t \bs\varepsilon_s') = \sigma_{ts}I_N\).
onde:
\(\bs X_t^*\): é a matriz de regressoras sem o intercepto;
\(\bs W_t^y\), \(\bs W_t^x\) e \(\bs W_t^u\): são as matrizes de ponderação, com dimensão \(n \times n\);
Modelo
Comumente, utiliza-se \(\bs W = \bs W_t^y = \bs W_t^x = \bs W_t^u\). Então, o modelo \(\eqref{eq: surs}\) pode ser reescrito como:
\[\bs A \bs y = \bs X \bs\beta + (I_{g} \otimes \bs W)\bs X^* \bs\theta + \bs u\] com \(\bs B \bs u = \bs\varepsilon\), \(\bs\varepsilon \sim \norm(\bs 0, \bs\Omega)\).
SUR-SDM: modelo espacial de Durbin, \(\lambda_t = 0\);
SUR-SDEM: modelo de erro espacial de Durbin; \(\rho_t = 0\);
SUR-SARAR: modelo de defasagem espacial com erros autorregressivos, \(\bs\theta_t = \bs0\).
Estimação
É necessário impor que a diagonal de \(\bs W\) seja preenchida por \(0\). Então, assumindo que os erros são normalmente distribuidos, basta maximizar a função log-verossimilhança, dada por:
\[\begin{split}
l(\bs y; \bs \eta) = - \frac{ng}{2} \ln(2\pi) - \frac{n}{2}\ln|\bs\Sigma| + g \sum_{t=1}^{g} \ln|\bs B_t| + g \sum_{t=1}^{g} \ln|\bs A_t| - \\-\frac{1}{2}(\bs A \bs y - \bs{\bar X}[\bs\beta; \bs\theta]')'\bs B' \bs\Omega^{-1} \bs B(\bs A \bs y - \bs{\bar X}[\bs\beta; \bs\theta]')
\end{split}\] onde \(\bs\eta'=[\bs\beta'; \bs\theta'; \rho_1;\dots;\rho_g; \lambda_1;\dots; \lambda_g; \sigma_{ij}]\), \(\bs\Lambda = diag(\lambda_1, \dots, \lambda_g)\), \(\bs\Gamma = diag(\rho_1,\dots,\rho_g)\), \(\bs{\bar X} = [\bs X; (T_g \otimes \bs W) \bs X^*]\), \(\bs A = I_{n g} - \bs\Gamma \otimes \bs W\), \(\bs B = I_{ng} - \bs\Lambda \otimes \bs W\) e \(\bs\Omega = \bs\Sigma \otimes I_n\).
Ao todo, temos que estimar \(2 \sum_{i=1}^g k_i + g + \frac{g(g+1)}{2}\) parâmetros 😨😱.
Aplicação 1
Conjunto de Dados
Banco de dados:
spc: este banco de dados ilustra um exemplo classico da curva de Phillips retirado de Anselin (\(1988\), pp. \(203\)–\(211\)).
Curva de Phillips
A curva de Phillips é uma teoria macroeconômica que mostra a relação entre a taxa de desemprego e a taxa de inflação no curto prazo. Ela foi criada pelo economista neozelandês Alban William Phillips, a partir de dados do Reino Unido entre \(1861\) e \(1957\).
spc
Objetivo do estudo: analisar como diferentes variáveis impactaram as mudanças nas taxas salariais de \(1981\) para \(1983\) em \(25\) condados do sudoeste de Ohio.
SMSA: variável dummy (\(1\) para condados metropolitanos, \(0\) caso contrário).
spc
Código
# Pacotes utilziados na aplicaçãorequire('gridExtra')require('ggplot2')require('dplyr')require('Rcpp')require('Gmisc')require('grid')require('spdep')require('sf')require('spsur')# Sintaxespcformula <- WAGE83 | WAGE81 ~ UN83 + NMR83 + SMSA | UN80 + NMR80 + SMSA# Criação do modelo SLMspcsur.slm <-spsurml(formula = spcformula, data = spc, type ='slm', method ='eigen', listw = Wspc)
neighbourhood matrix eigenvalues
Computing eigenvalues ...
Initial point: log_lik: 113.198 rhos: -0.472 -0.446
Iteration: 1 log_lik: 114.088 rhos: -0.506 -0.482
Iteration: 2 log_lik: 114.098 rhos: -0.506 -0.482
Iteration: 3 log_lik: 114.099 rhos: -0.505 -0.482
Time to fit the model: 0.72 seconds
Time to compute covariances: 0.21 seconds
type: slm (Spatial Lag Model), que considera \(\lambda_t = 0\) e \(\bs\theta_t= 0\), \(\forall t\).
method: “eigen”, utilizado para o cálculo do \(|A|\) e \(|B|\), sendo mais eficiente para amostras pequenas, mas computacionalmente caro para amostras grandes. Alternativa: realizar Fatoração LU ou Cholesky.
listw: lista de pesos espaciais, usada para definir vizinhanças em modelos espaciais.
spc
\(\rho\) é o parâmetro de autoregressão espacial, que indica a intensidade da dependência espacial nas variáveis dependentes. No nosso caso, \(\rho_1 = -0.505 < 0\) e \(\rho_2 = -0.482 < 0\) refletem uma dependência espacial negativa.
Bastou três iterações para a convergência do algoritmo, visto que não houveram mudanças significativas na função de verossimilhança e nos \(\rho\).
O teste de Breusch-Pagan (da matriz de covariância diagonal) é o responsável por testar se a diagonal da matriz de variância contemporânea é igual à \(0\). Note que \(\text{p-valor} = 0.0105 < 0.05\), logo, as equações possuem variâncias nulas, notável na saída da matriz de variância-covariância dos resíduos.
spc
A saída LMM é o multiplicador Lagrangiano marginal, usado para testar se foram omitidos efeitos espaciais decorrentes dos modelos estimados. Nesse caso, como \(\text{p-valor} = 0.477 > 0.05\) não temos evidências suficientes para rejeitar a hipótese nula em favor da alternativa, isto é, não há evidência de autocorrelação espacial dos resíduos.
Intervalos de confiança para os betas da primeira equação
Código
plot(spcsur.slm)
Intervalos de confiança para os betas da segunda equação
Intervalos de confiança para \(\rho\)
spc
Código
# Razão de verossimilhançasspcsur.slm <-spsurml(formula = spcformula, data = spc, type ='slm',listw = Wspc, control =list(trace =FALSE))spcsur.sdm <-spsurml(formula = spcformula, data = spc, type ='sdm',listw = Wspc, control =list(trace =FALSE))anova(spcsur.slm, spcsur.sdm)
É possível comparar se modelos de diferentes tipos são equivalentes com a função anova. Foi considerado o modelo SUR-SDM, que além de incluir a defasagem espacial das variáveis respostas, inclui para as preditoras.
Como \(\text{p-valor} = 0.4973 > 0.05\) não temos evidências suficientes para rejeitar a hipótese nula de que os modelos diferem.
Além disso, os critérios AIC e BIC penalizam modelos com mais parâmetros, busca-se valores menores desses critérios, indicando modelos mais parcimoniosos, isto é, uma explicação boa dos dados com o menor número possível de parâmetros.
spc
Código
# Teste de Wald, nesse caso estamos comparando se os interceptos são os mesmosR1 <-matrix(c(1, 0, 0, 0, -1, 0, 0, 0), nrow =1)b1 <-matrix(0, ncol =1)wald_betas(spcsur.slm, R = R1, b = b1)
Wald test on beta parameters
data: spc
Wald test = 0.39767, df = 1, p-value = 0.5283
Os interceptos não são estatisticamente diferentes entre si. Logo, uma especificação correta seria incluir o mesmo intercepto nas duas equações.
spc
Código
# adicionou um intercepto e fez um modelo restritoR1 <-matrix(c(1, 0, 0, 0, -1, 0, 0, 0), nrow =1)b1 <-matrix(0, ncol =1)spcsur.slm.restricted <-spsurml(formula = spcformula, data = spc, type ='slm', listw = Wspc,R = R1, b = b1, control =list(trace =FALSE))print(spcsur.slm.restricted)
coeff_1 pval_1 coeff_2 pval_2
(Intercept) 1.5786 0.000 1.5786 0.000
UN83 0.7946 0.003 NA NA
NMR83 -0.4873 0.064 NA NA
SMSA -0.0056 0.632 0.0029 0.906
UN80 NA NA -0.6193 0.120
NMR80 NA NA 0.7083 0.079
rho -0.5856 0.006 -0.3705 0.052
Código
# Teste que compara ro1 de ro2 DO MODELO RESTRITOR2 <-matrix(c(1, -1), nrow =1)b2 <-matrix(0, ncol =1)wald_deltas(spcsur.slm.restricted, R = R2, b = b2)
Wald test on spatial delta parameters
data: spc
Wald test = 16.932, df = 1, p-value = 3.874e-05
Confirmação de que a especificação com dois parâmetros diferentes da dependencia espacial estava correta.
Verificando normalidade dos resíduos
Código
residuos_eq1 <-residuals(spcsur.slm)[[1]]qqnorm(residuos_eq1, main ='QQ Plot dos Resíduos - Equação 1', col ='blue', pch =19, ylab ='Quantis Amostrais', xlab ='Quantis Teóricos')qqline(residuos_eq1, col ='red', lwd =2)
Verificando normalidade dos resíduos
Código
residuos_eq2 <-residuals(spcsur.slm)[[2]]qqnorm(residuos_eq2, main ='QQ Plot dos Resíduos - Equação 2', col ='blue', pch =19, ylab ='Quantis Amostrais', xlab ='Quantis Teóricos')qqline(residuos_eq2, col ='red', lwd =2)
Aplicação 2
NCOVR
O segundo conjunto de dados, retirado de Ballter et al. (\(2001\)), refere-se as taxas de homicídios em \(3085\) condados dos EUA para quatro anos (\(1960\), \(1970\), \(1980\) e \(1990\)).
Mínguez, R., López, F.A., & Mur, J. (2022). spsur: An R Package for Dealing with Spatial Seemingly Unrelated Regression Models. Journal of Statistical Software, 104(11), 1–43. doi: 10.18637/jss.v104.i11.
Anselin L (1988). Spatial Econometrics: Methods and Models. Studies in Operational Regional Science. Springer-Verlag, Dordrecht. doi: 10.1007/978-94-015-7799-1.
Baller, R.D. et al. (2001). Structural Covariates of US County Homicide Rates: Incorporating Spatial Effects. Criminology, 39(3), 561–588. doi: 10.1111/j.1745-9125.2001.tb00933.x.